R basic operations

# load library
suppressPackageStartupMessages({
  library(leaflet)
  library(sf)
  library(dplyr)
  library(lwgeom)
  library(purrr)
  library(DT)
  library(ggplot2)
  library(mapview)
  library(tmap)
  library(ggspatial)
})

## Warning in fun(libname, pkgname): GEOS versions differ: lwgeom has 3.12.1 sf
## has 3.12.2

## Warning in fun(libname, pkgname): PROJ versions differ: lwgeom has 9.3.1 sf has
## 9.4.1

Data

Data Source:
- NYC Borough Boundaries
- NYC Park Properties
- NYC Street Centerline

target_crs <- 2263

# load data
borough_boundaries <- st_read("D:/GIS/R_compare/Borough Boundaries/geo_export_43cec398-2125-4dd7-bd74-d3b7e1dfe9f5.shp")

## Reading layer `geo_export_43cec398-2125-4dd7-bd74-d3b7e1dfe9f5' from data source `D:\GIS\R_compare\Borough Boundaries\geo_export_43cec398-2125-4dd7-bd74-d3b7e1dfe9f5.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25559 ymin: 40.49613 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS:  WGS 84

parks <- st_read("D:/GIS/R_compare/Parks Properties_20241107/geo_export_d6f9bffe-7dce-4a9f-aae6-4ee8c4d82238.shp")

## Reading layer `geo_export_d6f9bffe-7dce-4a9f-aae6-4ee8c4d82238' from data source `D:\GIS\R_compare\Parks Properties_20241107\geo_export_d6f9bffe-7dce-4a9f-aae6-4ee8c4d82238.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2051 features and 33 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25609 ymin: 40.49449 xmax: -73.70905 ymax: 40.91133
## Geodetic CRS:  WGS 84

streets <- st_read("D:/GIS/R_compare/NYC Street Centerline (CSCL)/geo_export_03186438-85ec-4a9d-9f79-f75e3ba6fa7e.shp")

## Reading layer `geo_export_03186438-85ec-4a9d-9f79-f75e3ba6fa7e' from data source `D:\GIS\R_compare\NYC Street Centerline (CSCL)\geo_export_03186438-85ec-4a9d-9f79-f75e3ba6fa7e.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 121949 features and 34 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -74.25496 ymin: 40.49788 xmax: -73.70002 ymax: 40.9151
## Geodetic CRS:  WGS 84

streets <- streets[!st_is_empty(streets), ]

# Check and fix invalid geometries in parks
if (!all(st_is_valid(parks))) {
  parks <- st_make_valid(parks)  # Fix invalid geometries
}

R Basic Operations

1. Create a Map

Dynamic maps: Packages such as leaflet and mapview in R support the creation of interactive maps that can be dynamically displayed on web pages. In ArcGIS Pro, similar interactivity requires the use of ArcGIS Online.
Flexible custom layers and layouts: ggplot2 and tmap packages allow users to freely customize map layouts and layouts, suitable for complex multi-layer drawing.

leaflet

Leaflet is a powerful open-source JavaScript library used to create interactive maps. The R package ‘leaflet’ provides a convenient interface to integrate Leaflet’s capabilities within R, allowing users to build dynamic maps with ease. The package supports various map elements like tiles, markers, polygons, and popups, and is designed to be simple, efficient, and user-friendly.

# R: Dynamic Maps

# Create a Leaflet map with a basemap, legend, title, and scale
leaflet(parks) %>%
  addTiles() %>%  # Add online basemap
  addPolygons(color = "green", weight = 1, fillOpacity = 0.5) %>%
  addLegend(position = "bottomright", colors = "green", labels = "NYC parks", title = "Legend") %>%
  addScaleBar(position = "bottomleft") %>%
  addControl("<strong>NYC parks Map</strong>", position = "topright")

Advantages:

• Interactive: The map generated by Leaflet can be panned and zoomed, supports click events and pop-up windows, and is very suitable for displaying interactive maps online.

• Lightweight and fast: Leaflet is a lightweight library that loads quickly and runs smoothly even on most browsers.
  

• Open source and free: Compared with the commercial license fee of ArcGIS Pro, Leaflet is completely free and has rich community support and plug-in library.

Disadvantages:

• Limited analysis function: Leaflet is mainly used to display maps and lacks spatial analysis functions. If complex spatial analysis is required, it is necessary to combine other tools or data preprocessing.
• Limited rendering quality: Leaflet’s rendering quality is relatively simple, which is suitable for basic map display, but not suitable for making high-quality and complex maps.It support only crs.
• Limited CRS: Because Leaflet uses longitude and latitude coordinates to express coordinates, the coordinate system used in Leaflet is WGS84, which is EPSG:4326. Therefore, whether calculating the row and column numbers of tiles or calculating the final position of tiles on the map, the corresponding coordinate system definition is required for conversion.

mapview

Mapview is a fast visualization package in R, specifically for displaying geospatial data, based on Leaflet.

# Create an interactive map using mapview
map <- mapview(parks, legend = TRUE, layer.name = "NYC Parks", map.types = c("OpenStreetMap", "Esri.WorldImagery"))

# Show the map
map

Advantages:

• Get started quickly: mapview can automatically load and display simple spatial data (such as points, lines, and surfaces) without complex code.
  
• Good interactivity: Inheriting the interactive features of Leaflet, users can pan, zoom, and click to view data attributes, which is suitable for quickly browsing spatial data.
  
• Support multiple layers: Multiple layers can be easily superimposed, which is more suitable for quickly viewing multiple data layers than Leaflet.

Disadvantages:

• Limited analysis functions: It also lacks complex spatial analysis functions and is mainly used for data display.
  
• Slightly less customizable: Although it can display multiple layers, it is not as good as ArcGIS Pro in terms of custom styles and detail adjustments.

ggplot2

ggplot2 is one of the most commonly used visualization libraries in R, suitable for making high-quality static maps.

# Create a static map using ggplot2
ggplot(borough_boundaries) +
  geom_sf(aes(fill = boro_name), color = "black") +  # 'boro_name' represents borough names
  scale_fill_discrete(name = "Borough") +  # Add legend for borough names
  labs(
    title = "NYC Borough Boundary",
    caption = "Data Source: NYC OpenData"
  ) +
  # Add scale bar and north arrow
  annotation_scale(location = "bl", width_hint = 0.3) +  # Add scale bar at the bottom left
  annotation_north_arrow(location = "tl", which_north = "true", style = north_arrow_fancy_orienteering) +  # Add north arrow at the top left
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    plot.caption = element_text(hjust = 0.5, size = 10),
    legend.position = "bottom"
  )

Advantages:

• High-quality output: ggplot2 supports rich customization and style control, and can create publishing-level static maps.
• Good data integration: ggplot2 can seamlessly integrate other data analysis tools in R, making it easy to visualize data after operation.
• Suitable for making complex statistical maps: such as stratification, classification, heat map, etc., and can flexibly use variables such as color, shape and transparency.

Disadvantages:

• Lack of interactivity: ggplot2 generates static maps and cannot be interactively operated.
• Geographic data processing is slightly cumbersome: Compared with ArcGIS Pro, other R packages (such as sf) are required to process spatial data, and the process may be more complicated.

tmap

tmap is a package in R specifically for map making, suitable for static and interactive map production.

# Enable automatic fixing of invalid polygons (if needed)
tmap_options(check.and.fix = TRUE)

# Set tmap to static mode
tmap_mode("plot")

## tmap mode set to plotting

# Create a static map with a basemap
tm_shape(borough_boundaries) + 
  tm_polygons(title = "Borough Boundaries", col = "lightgrey", border.col = "black", legend.show = TRUE) +
  tm_tiles("OpenStreetMap") +  # Add a static basemap after all layers
  tm_layout(
    title = "NYC Map with Boroughs, Parks, and Streets",
    title.size = 1.5,
    title.position = c("center", "top"),
    legend.position = c("left", "bottom"),
    legend.bg.color = "white",
    legend.frame = TRUE,
    frame = FALSE,
    bg.color = "lightblue",
    inner.margins = c(0.05, 0.05, 0.05, 0.05)
  ) +
  tm_scale_bar(position = c("left", "bottom")) +
  tm_compass(type = "8star", position = c("right", "top")) +
  tm_credits("Data Source: NYC Open Data", position = c("right", "bottom"))

# Set tmap to interactive mode
tmap_mode("view")

## tmap mode set to interactive viewing

# Create an interactive map
tm_shape(borough_boundaries) + 
  tm_polygons(title = "Borough Boundaries", col = "lightgrey", border.col = "black") +
  tm_layout(
    title = "NYC Map with Boroughs, Parks, and Streets",
    title.size = 1.5,
    title.position = c("center", "top"),
    legend.position = c("left", "bottom"),
    legend.bg.color = "white",
    legend.frame = TRUE,
    frame = FALSE,
    bg.color = "lightblue",
    inner.margins = c(0.05, 0.05, 0.05, 0.05)
  ) +
  tm_basemap("OpenStreetMap") +  # Add an online basemap
  tm_scale_bar(position = c("left", "bottom"))

## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.

Advantages:

• Mode switching: tmap supports static and interactive modes (“plot” and “view” modes), which can balance high-quality map output and interactivity.
  
• Professional mapping function: tmap supports a variety of mapping options and is suitable for making thematic maps, such as zoning maps, bubble maps, etc.
  
• Easy to customize: It provides a variety of customization options, close to the mapping quality of ArcGIS, and can control details such as color, annotation, legend position, etc.

Disadvantages:

• Not suitable for complex analysis: Although tmap performs well in mapping, its spatial analysis capabilities are not as good as ArcGIS Pro.
  
• Slightly slower: The rendering speed is slow for large data sets or complex layers, which is suitable for static map generation rather than real-time display.

2. Extract Toolset - Clip

Multi-layer clipping operations: R supports chaining operations, which allows you to clip multiple layers in sequence in a single code execution. In ArcGIS Pro, this requires multiple separate operations to complete.

Get the parks within Brooklyn.

# 1. Get the Brooklyn boundary
brooklyn_boundary <- borough_boundaries[borough_boundaries$boro_name == "Brooklyn", ]
brooklyn_boundary <- brooklyn_boundary[st_is_valid(brooklyn_boundary), ]
brooklyn_street<-streets[streets$borocode == 3, ]
brooklyn_parks<-parks[parks$borough == "B", ]

brooklyn_parks <- st_make_valid(brooklyn_parks)
brooklyn_street <- st_make_valid(brooklyn_street)

# 3. Perform multi-layer clipping to get parks near streets

parks_in_Brooklyn <- st_intersection(parks, brooklyn_boundary)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

map1 <- mapview(brooklyn_parks, color = "blue", fill = "blue", alpha.regions = 1, lwd = 2, layer.name = "Brooklyn Parks by selection") +
  mapview(parks_in_Brooklyn, color = "green", fill = "green", alpha.regions = 1, lwd = 2, layer.name = "Brooklyn Parks by clipping") +
  mapview(brooklyn_boundary, color = "red", fill = NA, lwd = 2, layer.name = "Brooklyn Boundary")


map1

3. Extract Toolset - Select

Flexible conditional expressions: R’s dplyr package can use more complex conditional expressions and functions to select data, supporting any custom conditions. In ArcGIS Pro, using the GUI to make such complex selections requires more steps and multiple conditional expressions.

Select areas in Brooklyn that have green space larger than a certain value.

brooklyn_boundary <- borough_boundaries %>% filter(boro_name == "Brooklyn")

# 2. parks in Brooklyn
suppressWarnings({parks_in_brooklyn <- st_intersection(parks, brooklyn_boundary)})

# 3. add area column
parks_in_brooklyn <- parks_in_brooklyn %>%
  mutate(area_m2 = st_area(geometry))

# 4. select parks whose area > 1.5millions square meters
large_parks <- parks_in_brooklyn %>% filter(area_m2 > units::set_units(1500000, "m^2"))

# Plot the map and add a custom legend
ggplot() +
  # Plot Brooklyn boundary (light grey)
  geom_sf(data = brooklyn_boundary, aes(fill = "Brooklyn Boundary"), color = "lightgrey", size = 0.5, show.legend = "fill") +
  
  # Plot large parks
  geom_sf(data = large_parks, aes(fill = "large parks"), color = "green", size = 0.1, show.legend = "fill") +
  
  # Add title and legend
  labs(
    title = "Parks in Brooklyn larger than 1.5 M (2024)",
    caption = "Data Source: NYC Open Data",
    fill = "Legend"
  ) +
  
  # Set color mapping
  scale_fill_manual(values = c(
    "Brooklyn Boundary" = "lightgrey",
    "large parks" = "green"
  )) +
  
  # Control border display in legend only
  guides(fill = guide_legend(override.aes = list(color = "white", size = 0.5))) +
  
  # Set theme style
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 12),
    plot.caption = element_text(hjust = 1, size = 10),
    legend.position = "bottom"
  )

4. Extract Toolset - Split

Complex split operations: R supports splitting based on arbitrary geometry, and the st_split() function can split a layer based on custom shapes or other conditions. In ArcGIS Pro, splitting by arbitrary shapes is not always straightforward and requires complex layer editing. Compared with ArcGIS Pro, we can create a loop to split items in a list.

Use street data to segment Brooklyn park area.

large_parks <- large_parks[!st_is_empty(large_parks), ]

split_parks <- map(large_parks$geometry, ~ st_split(.x, streets))

split_parks_geom <- map(split_parks, ~ st_collection_extract(.x, "POLYGON"))

split_parks_combined <- do.call(c, split_parks_geom) %>%
  st_sf()

map2 <- mapview(brooklyn_boundary, 
               col.regions = "lightgrey", 
               color = "darkgrey", 
               layer.name = "Brooklyn Boundary") +
       mapview(split_parks_combined, 
               col.regions = "green", 
               color = "darkgreen", 
               layer.name = "Split Parks Combined")

## Warning in cbind(`Feature ID` = fid, mat): number of rows of result is not a
## multiple of vector length (arg 1)

map2

5. Spatial Join

Join long streets’ attribute to the 2 largest park in Brooklyn

long_streets <- streets %>%
  filter(st_length(geometry) > units::set_units(500, "m"))
joined_parks_streets <- st_join(large_parks, long_streets, join = st_intersects)

datatable(
  joined_parks_streets,
  options = list(pageLength = 5, autoWidth = TRUE, scrollX = TRUE),
  caption = "Joined Parks and Streets Attributes Table"
)

6. Buffer

Dynamic buffer: R supports the creation of dynamic buffers based on attributes, such as adjusting the buffer distance based on the attribute value of an object. In ArcGIS Pro, this attribute-based dynamic buffer needs to be implemented through a combination of attribute field settings and the buffer tool.

Find the buffer for the long streets based on their width.

buffered_polygon <- long_streets %>% st_buffer(dist = long_streets$st_width)

# Plot the map and add a custom legend
ggplot() +
  # Plot Brooklyn boundary (light grey)
  geom_sf(data = borough_boundaries, aes(fill = "NYC Boundary"), color = "darkgrey", size = 0.5, show.legend = "fill") +
  
  # Plot buffer
  geom_sf(data = buffered_polygon, aes(fill = "buffer zone"), color = "green", size = 0.1, show.legend = "fill") +
  
  # Add title and legend
  labs(
    title = "Buffer for Long Street in NYC",
    caption = "Data Source: NYC Open Data",
    fill = "Legend"
  ) +
  
  # Set color mapping for fill and border
  scale_fill_manual(values = c(
    "NYC Boundary" = "lightgrey",
    "buffer zone" = "green"
  )) +
  
  guides(fill = guide_legend(override.aes = list(color = "darkgrey", size = 0.3))) +
  
  # Set theme style
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    plot.caption = element_text(hjust = 1, size = 10),
    legend.position = "bottom"
  )